Simple Option Pricing

Introduction

Introduction

So far you will be familiar with the binomial model and the Black Scholes model for option pricing

Fundamentally as a result of the Cameron Martin Girsanov Theorem and the Central Limit Theorem both these methods are effectively the same thing

As the number of steps tends to $\infty$ the binomial model tends to the Black Scholes model

The Monte Carlo methods are also fundamentally the same thing - but with the Monte Carlo methods we can start from a binomial perspective or from a Black Scholes perspective and gain insight into how all the methods are effectively equivalent

Additionally with the Monte Carlo method it is much easier to extend the method to deal with more advanced modelling of equity prices with features such as skew and stochastic volatility

$\mathbb{P}$ and $\mathbb{Q}$ Review

$\mathbb{P}$ and $\mathbb{Q}$ Review

Remember when you first looked at the binomial model

You started with an uptick and a downtick

Then do some algebra and the formula that emerges is $f_{n} = e^{-r\Delta t} \times E_\mathbb{Q} \left(f_{n+1}\right)$, that is effectively an expected value over some 'other' probability measure '$\mathbb{Q}$'

Review of CRR model

We start with a very simple model of the price evolution of our underlying stock

The stock has price $S_0$ at time 0 and the model advances in units of time $\Delta t$

The term of the option is $T$

We assume the volatility of the stock is $\sigma$, which leads to a natural choice of $u$ and $d$ of:

$u=e^{\sigma \sqrt{\Delta t}}$ and $d=\frac{1}{u}$

and $q=\frac{e^{r \Delta t}-d}{u-d}$ - from arbitrage free construction

The value at each node is then worked out from the subsequent nodes with the formula: $f_{n} = e^{-r\Delta t} \times E_\mathbb{Q} \left(f_{n+1}\right)$

These formulas are derived by considering a hedging portfolio of stocks and cash

The important thing about $u$ and $d$ is the gap between them as the Girsanov Theorem can sort out the drift but not the volatility

The $\sqrt{\Delta t}$ is to make the unit time volatility $\sigma$

So simply by changing the probabilities of the uptick and the downtick we can express the market price of an option (or any other financial asset) as an expected value

This extends through the Cameron Martin Girsanov Theorem to the full Black Scholes Model.

An illustration of this can be seen in this graphical simulation

Binomial Monte Carlo

Binomial Monte Carlo

Weighted Probabilities

Much of the mathematics of the binomial method has to be reproduced to use the Monte Carlo method

although Monte Carlo methods are often simpler than other methods we CANNOT

All we can do in the Monte Carlo method is simplify the algorithm by being able to pick the up tick and the down tick randomly rather than having to write the code to go through every single permutation

For given values of $S_0$, $K$, $T$, $\Delta T$ and $\sigma$ can you write a Monte Carlo version of the binomial model

Can you explain how this is simpler to code than the normal binomial model

There is just one path for each iteration

What other advantages does it have?

You do not run the risk of running out of memory especially for those options which are not recombining

Distribution Shift

Fundamentally the CMG is about $\frac{d \mathbb{Q}}{d \mathbb{P}}$ but the visible impact of this is to appear to shift the distribution from the expected return on equity to the risk free rate as we have seen in the above demonstration

Can you change your model to reflect this directly in the VBA

Hint: Preserve the ratio of $u:d$ and solve $q=0.5$

Solution:

$u=C \times e^{\sigma \sqrt{\Delta t}}$ and $d=C \times e^{-\sigma \sqrt{\Delta t}}$

$q=\frac{e^{r \Delta t}-d}{u-d} = 0.5$

$q=\frac{e^{r \Delta t}-C \times e^{-\sigma \sqrt{\Delta t}}}{C \times e^{\sigma \sqrt{\Delta t}}-C \times e^{-\sigma \sqrt{\Delta t}}} = 0.5$

$C=\frac{e^{r \Delta t}}{e^{-\sigma \sqrt{\Delta t}} + 0.5 \times e^{\sigma \sqrt{\Delta t}} - 0.5 \times e^{-\sigma \sqrt{\Delta t}}}$

$C=\frac{2 \times e^{r \Delta t}}{e^{\sigma \sqrt{\Delta t}} +e^{-\sigma \sqrt{\Delta t}}}$

Multiple Step Lognormal

Multiple Step Lognormal

It is a necessary limitation of the binomial model to only have discrete state spaces.

This clearly means an uptick and a downtick in the case of the binomial model but a basic extension is the trinomial model which has three ticks from each state

However for a Monte Carlo approach we are not limited to discrete state spaces.

Why

because we do not store each state

Therefore it is possible (and sensible) to extend the model by allowing the return at each point to be generated fully from a lognormal distribution

Extend your model to allow for full lognormal steps

Monte Carlo Options

How many steps is it now sensible to have to get to an accurate result (You may assume sufficient Monte Carlo runs).

ONE as we only needed multiple steps under the binomial model to allow time for the central limit theorem to get us to a lognormal distribution

For simple options with lognormal steps ONE step is the Monte Carlo equivalent of the Black Scholes model and will converge to the Black Scholes result as the number of runs becomes large

One Step Monte Carlo

One Step Monte Carlo

Test your model using just one step

Can you reconcile the mathematics of what we have done to the proof of the Black Scholes formula

The Monte Carlo now is simply the final stage of the Black Scholes proof - that is the integration at the end once we have established that the arbitrage free value is genuinely the expected value over $\mathbb{Q}$

We still need Ito's Lemma

Ito's Lemma is Still Necessary

If you have followed the steps above you will probably expect to have replicated Black Scholes results with a one-step lognormal model.

However (unless you have thought ahead to this section) it is highly likely you will have made a mistake and your answers will be too high

The mistake we have made is to assume that the solution of the integral equation: $\frac{dS_t}{S_t}=r dt + \sigma dW_t$ is $S_t=S_0 e^{rt+\sigma W_t}$

When in fact it is $S_t=S_0 e^{(r-0.5 \sigma^2)t+\sigma W_t}$

Proof

$\frac{dS_t}{S_t}=r dt + \sigma dW_t \implies dS_t=r S_t dt + \sigma S_t dW_t$

Now look at Ito's Lemma:

$df(t, X_t) = \left[\frac{df}{dt} + \mu (t, X_t) \frac{df}{dx} + \frac{1}{2} \sigma^2 (t, X_t) \frac{d^2f}{dx^2} \right] dt + \sigma(t, X_t) \frac{df}{dx} dW_t$

Where $dX_t=\mu(t,X_t)dt + \sigma(t, X_t) dW_t$

We set this up as follows:

Rewrite with $S_t$ instead of $X_t$

Then our $\mu$ is $r S_t$ and our $\sigma$ is $\sigma S_t$

The we use the function $f(t, S_t)=ln(S_t)$ (there is no way of knowing this other than it works)

so $df(t, S_t) = \left[\frac{df}{dt} + \mu (t, S_t) \frac{df}{ds} + \frac{1}{2} \sigma^2 (t, S_t) \frac{d^2f}{ds^2} \right] dt + \sigma(t, S_t) \frac{df}{ds} dW_t$

then $df(t, S_t) = \left[\frac{df}{dt} +r S_t \frac{df}{ds} + \frac{1}{2} \sigma^2 S_t^2 \frac{d^2f}{ds^2} \right] dt + \sigma S_t \frac{df}{ds} dW_t$

Now we need to calculate our derivatives:

$f(t, S_t)=ln(S_t) \implies \frac{df}{dt}=0, \frac{df}{ds}=\frac{1}{S_t}, \frac{d^2f}{ds^2}=\frac{-1}{S_t^2}$

So putting it all together:

$d ln(S_t)= \left[0 +r S_t \frac{1}{S_t} + \frac{1}{2} \sigma^2 S_t^2 \frac{-1}{S_t^2} \right] dt + \sigma S_t \frac{1}{S_t} dW_t$

Simplify

$d ln(S_t)= \left[r - \frac{1}{2} \sigma^2 \right] dt + \sigma dW_t$

Integrate:

$\displaystyle\int d ln(S_t)= \displaystyle\int (r - \frac{1}{2} \sigma^2) dt +\displaystyle\int \sigma dW_t$

$ln(S_t)= (r - \frac{1}{2} \sigma^2)t + \sigma W_t + ln S_0$

$\therefore S_t=S_0 e^{(r-\frac{\sigma^2}{2})t + \sigma W_t}$

Code for Different Methods

Code for Different Methods

We can now see below the code for the three different methods

Option Explicit
Function max(x As Double, y As Double) As Double
  If x > y Then
    max = x
   Else
    max = y
   End If
End Function

Function CRR(S0, K, r, T, sigma As Double, n As Integer) As Double
  Dim S() As Double
  ReDim S(0 To n, 0 To 2 ^ n - 1)
  Dim f() As Double
  ReDim f(0 To n, 0 To 2 ^ n - 1)
  Dim delta_t, u, d As Double
  Dim j As Double
  Dim q As Double
  Dim node_x, node_y As Integer
  delta_t = T / n
  u = Exp(sigma * delta_t ^ 0.5)
  d = Exp(-sigma * delta_t ^ 0.5)
  q = (Exp(r * delta_t) - d) / (u - d)
  S(0, 0) = S0
  For node_x = 1 To n
    For node_y = 0 To 2 ^ node_x - 1
      If node_y Mod 2 = 0 Then
        S(node_x, node_y) = S(node_x - 1, node_y / 2) * u
       Else
        S(node_x, node_y) = S(node_x - 1, (node_y - 1) / 2) * d
       End If
     Next node_y
   Next node_x
   
  For node_y = 0 To 2 ^ n - 1
    f(n, node_y) = max(S(n, node_y) - K, 0)
   Next node_y
   
 For node_x = n - 1 To 0 Step -1
   For node_y = 0 To 2 ^ node_x - 1
     f(node_x, node_y) = Exp(-r * delta_t) * (q * f(node_x + 1, node_y * 2) + (1 - q) * f(node_x + 1, node_y * 2 + 1))
    Next node_y
  Next node_x
CRR = f(0, 0)
End Function

Function Monte_Carlo_Binomial(S0, K, r, T, sigma As Double, n, runs As Integer) As Double
  Dim delta_t, u, d As Double
  Dim j As Long
  Dim q As Double
  Dim S As Double
  Dim sum As Double
  Dim payoff As Double
  Dim run As Long
  delta_t = T / n
  u = Exp(sigma * delta_t ^ 0.5)
  d = Exp(-sigma * delta_t ^ 0.5)
  q = (Exp(r * delta_t) - d) / (u - d)
  sum = 0
  For run = 1 To runs
    S = S0
    For j = 1 To n
      If (Rnd() < q) Then
        S = S * u
       Else
        S = S * d
       End If
     Next j
    payoff = max(S - K, 0)
    sum = sum + payoff
   Next run
  Monte_Carlo_Binomial = Exp(-r * T) * sum / runs
End Function

Function Monte_Carlo_LogNormal(S0, K, r, T, sigma As Double, n, runs As Integer) As Double
  Dim delta_t, u, d As Double
  Dim j As Long
  Dim q As Double
  Dim S As Double
  Dim jump As Double
  Dim sum As Double
  Dim payoff As Double
  Dim run As Long
  delta_t = T / n
  sum = 0
  For run = 1 To runs
    S = S0
    For j = 1 To n
      jump = Exp(Application.NormSInv(Rnd()) * sigma * Sqr(delta_t) + (r-0.5*sigma^2) * delta_t)
      S = S * jump
     Next j
    payoff = max(S - K, 0)
    sum = sum + payoff
   Next run
  Monte_Carlo_LogNormal = Exp(-r * T) * sum / runs
End Function